%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Appendix Figure 8: Empirical derivative estimates in low-income neighborhoods out to 600m
%%
%% NOTE: Run rent_processing_macros.do before running this .m file
%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
addpath('CODE DIRECTORY')
addpath('DATA DIRECTORY')
addpath('OUTPUT DIRECTORY')

% Use 600m file %
cd 'DATA DIRECTORY'
file = dir('listing_data_w600m.csv');

cd 'OUTPUT DIRECTORY'

name_suffix=extractBetween(file.name,14,18);
perimeter=extractBetween(name_suffix,2,4);
perimeter=str2num(perimeter{1})/1000;

% Pull in listing data 
m=readtable(file.name,'ReadVariableNames',true);

%% Load main dataset and impose sample restrictions %% 
subsample=m(m.listing_bedrooms>0 & m.listing_bedrooms<4 & m.listing_bathrooms~=0 & ((m.building_built>=2015 & m.building_built<=2016)) & m.income_g==1,:);

%%%%%%%%%%%%%%%%%%%%%%
% Declare parameters %
%%%%%%%%%%%%%%%%%%%%%%

kappa=5;
theta_angle_n=0.45;
mi_band=0.2;
t_band=3;

%Count new building sites 
[rows,cols]=size(subsample);
subsample.building_id=categorical(subsample.building_id);
building_id_list=unique(subsample.building_id);
[total_sites,total_cols]=size(building_id_list);
data_points=zeros(total_sites,1);

%Initialize Variables
fips=zeros(total_sites,2);
smooth_out = cell(total_sites,1);

%Loop over building sites
for i=1:total_sites

    %%% Use the building_id_list vector as an index for each new building

    %Pull out listing price data around specific building site
    %listing prices
    id=char(building_id_list(i));
    subdata=subsample(subsample.building_id==id,:);

    % Generate residuals to partial out bed and baths %
    [rows,cols]=size(subdata);
    if rows>30
        bed1=[subdata.listing_bedrooms==1];
        bed2=[subdata.listing_bedrooms==2];
        bed3=[subdata.listing_bedrooms==3];
        bath1=[subdata.listing_bathrooms==1];
        bath2=[subdata.listing_bathrooms==2];
        bath3=[subdata.listing_bathrooms>=3];
        halfbath=[mod(subdata.listing_bathrooms,1)==0.5];
        X=[ones(rows,1) bed1 bed2 bed3 bath1 bath2 bath3 halfbath];

        b = regress(log(subdata.rent), X);
        Yhat = X*b;
        price = log(subdata.rent) - Yhat;

        year_built=subdata.building_built(1);

        %Year building built relative to 2011 (year before sample
        %starts)
        yr_all=subdata.building_built-2011;

        %Distance to new building site
        dist2=subdata.dist/1000;

        %Lat-lon of apartment listing
        coord=[subdata.listing_lat subdata.listing_lon];

        %Years of rent relative to new building
        %construction
        t=(subdata.listing_year-2011)-yr_all;

        %Neighborhood type
        data_points(i)=size(subdata,1);

        %Estimate empirical derviatives
        d=zeros(size(t,1),1);
        for ii=1:size(price)
            d(ii)=first_derivative_r_prepost(coord(ii,:),dist2(ii),t(ii),price,dist2,t,coord,theta_angle_n,kappa);
        end

        dd_clean=d(d~=0);
        coord_clean=coord(d~=0,:);
        dist_clean=dist2(d~=0);
        t_clean=t(d~=0);

        if ~isempty(t_clean>0) 
        %bandwidth choice        
                h=[mi_band t_band];

                phat=zeros(10,7);
                pphat=zeros(1,7);
                %NW Kernal smooth empirical derivatives
                for ii=1:10
                    for j=-3:3

                        x=[dist_clean t_clean];
                        xeval=[ii*.1*perimeter j+.01];

                        step1a=repmat(xeval,size(x,1),1);
                        step1b=repmat(h,size(x,1),1);
                        step1c=(x-step1a);
                        u=step1c./step1b;

                        k=prod(.75*(1-(u.^2)).*(abs(u)<=1),2);
                        k=k.*((x(:,2).*xeval(2))>0);

                        dd_out=[dd_clean'*k/sum(k); NaN ; NaN]; 

                        pphat(j+4)=dd_out(1);      
                    end
                    %Difference relative to year t-1 (pre-post)
                    avg=pphat(3);
                    pphat=pphat-avg;
                    phat(ii,:)=pphat;

                end
                %Final estimate of for this specific building site
                smooth_out{i}=phat;
            end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% Generate the Underlying Data for the Mesh Plot  %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

est_all=[];
bad=[];
save_pics=1;

points=data_points;
est=smooth_out;

est_all=[est_all; est];


for jk=1:size(est_all,1)
    est=est_all{jk};
    if size(est,1)>0  && (sum(sum(isnan(est(1:10,2:7))))<15)          
        est_all{jk}=est;
    else
        bad=[bad; jk];
    end
end
    est_all(bad)=[];

ct=zeros(10,7);
eest=zeros(10,7);

%actual estimate
  for i=1:size(est_all,1)
        data=est_all{(i)};   
        points=1;
        wt=points*(~isnan(data));
        data(isnan(data))=0;
        ct=ct+wt;
        eest=eest+(points*data);
  end

est_real=eest./ct;
ct_real=ct;

est_real=.1*est_real;

 %Derivative Graph
for i=1:7
    for j=1:10        
        %real estimates
        pp_real(j,i)=10*((est_real(j,i)));
        pp_lev_real(j,i)=-sum(sum(est_real(j:end,i)));
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Print the Mesh Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fig_save_title=['appendix_figure_8.png'];

base_mesh_value=perimeter/10;
%Appendix Figure 8
if save_pics==1
    figure
    mesh(-3:3,base_mesh_value:base_mesh_value:perimeter,pp_lev_real)
    xlabel('Years Since Construction')
    ylabel('Kilometers to Building Site')
    zlabel('Log Rent')
    ylim([0,perimeter])
    saveas(gcf,fig_save_title)
end


